## 

rm(list = ls())

## 

library(rdrobust)
library(tidyverse)
library(pbapply)

## Declare covars

covs <- c('turnout_party_2009_btw', 'soz_vers_beschaeftigte_share',
          'pop_density_km2', 'migration_out_share')

## Get federal election data

bt <- readRDS("data/data_federal.rds") %>%
  filter(!is.na(treated))

## List of outcomes

outcomes <- c('turnout_party',
              'agg_left_party', 
              'agg_center_party', 
              'agg_right_party',
              "right_total_party", 
              "left_total_party",
              "cdu_csu_party", 
              "spd_party", 
              "greens_party", 
              "fdp_party",
              "left_party", 
              "afd_party", 
              "other_party",
              "incumbent_party")

## Source helper function to tidy rdrobust output

source("code/tidy_rd.R")

## First differences

diff_df <- pblapply(outcomes, function(o) {
  out <- bt %>%
    filter(year > 2012) %>%
    pivot_wider(values_from = o, names_from = 'year', id_cols = 'ags',
                names_prefix = 'o') %>%
    mutate(diff = o2017  - o2013) %>%  dplyr::select(ags, diff) 
  ## Rename
  colnames(out)[2] <- o
  
  ## Return this
  out
}) %>%
  reduce(left_join) %>%
  left_join(bt %>% dplyr::select(ags, pop_dec_09, applies_census,
                                 one_of(covs), state_id) %>%
              distinct(ags, .keep_all = T)) %>%
  mutate(runvar = (pop_dec_09 * -1) + 10000) %>%
  filter(applies_census == 1)

## Declare list of bandwidths

bw_vec <- c(seq(1000, 5000, 1000)) %>% sort(decreasing = F)

## Estimate (this estimates effects for different BWs for all outcomes)
## This is mainly used for bandwidth sensitivity checks 

out_rd <- pblapply(outcomes, function(o) {
  
  # This iterates over bandwidths
  
  lapply(bw_vec, function(h) {
    
    # This iterates over whether or not to exclude B-W
    
    lapply(c(T,F), function(exclude_state) {
      
      ## Declare exclusion of B-W (or not)
      
      if (exclude_state) {
        subset_select <- !diff_df$state_id == '08' 
      } else {
        subset_select <- rep(T, nrow(diff_df))
      }
      
      ## Non-missing obs.
      
      complete_obs <- complete.cases(diff_df[, c(covs, o, 'runvar')]) &
        subset_select
      
      ## Estimate 
      
      out <- rdrobust(y = diff_df[complete_obs, o] %>% pull(!!o), 
                      x = diff_df[complete_obs, 'runvar'] %>% pull(runvar), 
                      c = 0, 
                      h = h, covs = diff_df[complete_obs, covs])
      
      ## Tidy and return
      
      out2 <- out %>% 
        tidy_rd(se_nr = 3) %>% 
        mutate(outcome = o, 
               exclude_state = ifelse(exclude_state, 
                                      'State excluded', 'All states')) %>% 
        dplyr::rename(se = std.error, pval = p.value, 
                      bw = bw_left_h,
                      bw_bias = bw_left_b) %>% 
        mutate(conf.low90 = estimate - qnorm(0.95)*se,
               conf.high90 = estimate + qnorm(0.95)*se)
      
      out2
    }) %>% reduce(rbind) 
  }) %>% reduce(rbind) 
}) %>% reduce(rbind)

## Next, estimate main effect w/ optimal BWs

out_rd_opt <- pblapply(outcomes, function(o) {
  
  ## Iterate over whether to exclude B-W
  
  lapply(c(T,F), function(exclude_state) {
    
    if (exclude_state) {
      subset_select <- !diff_df$state_id == '08' 
    } else {
      subset_select <- rep(T, nrow(diff_df))
    }
    
    ## Estimate
    
    out <- rdrobust(y = diff_df[subset_select, o] %>% pull(!!o), 
                    x = diff_df$runvar[subset_select], c = 0,
                    covs = diff_df[subset_select, covs])
    
    ## Tidy and return
    
    out2 <- out %>% 
      tidy_rd(se_nr = 3) %>% 
      mutate(outcome = o,
             exclude_state = ifelse(exclude_state, 
                                    'State excluded', 'All states')) %>% 
      dplyr::rename(se = std.error, pval = p.value, 
                    bw = bw_left_h,
                    bw_bias = bw_left_b) %>% 
      mutate(conf.low90 = estimate - qnorm(0.95)*se,
             conf.high90 = estimate + qnorm(0.95)*se)
    
    out2
  }) %>% reduce(rbind)
}) %>%
  reduce(rbind)

## Save results

write_rds(out_rd_opt, 'data/rdd_federal_results.rds')
write_rds(out_rd, 'data/rdd_federal_results_sensitivity.rds')
